#include<iostream>
using namespace std;
#include<fstream>
#include<cmath>

#include"library.h"

int main(){

	int i,j,k,n=4;
	double y;
	double A[4][4] = { {1,1,1,1}, {1,2,2,2}, {1,2,3,3}, {1,2,3,4} };

	double **a = new double *[n];
	for(i=0;i<n;i++){
		a[i] = new double [n];
		for(j=0;j<n;j++){
			a[i][j] = A[i][j];
		}
	}

	double **a_reverse = new double *[n];
	for(i=0;i<n;i++) a_reverse[i] = new double [n];

	double *b = new double [n];
	int *pivot = new int [n];
	int sign = 1;

	ludcmp(4, a, pivot, sign);

	for(j=0;j<n;j++){ //calculate the jth column of a_reverse	
		for(i=0;i<n;i++){
			if(i==j) b[i]=1;
			else b[i] = 0;
		}
		ludcmp_solve(n, a, pivot, b);
		for(i=0;i<n;i++){
			a_reverse[i][j] = b[i];
		}
	}
	cout<<"a_reverse:"<<endl;
	for(i=0;i<n;i++){
		for(j=0;j<n;j++){
			cout<<a_reverse[i][j]<<"\t";
		}
		cout<<endl;
	}
	cout<<"test if a * a_reverse = E"<<endl;
	for(i=0;i<n;i++){
		for(j=0;j<n;j++){
			y=0;
			for(k=0;k<n;k++){
				y += A[i][k] * a_reverse[k][j];
			}
			cout<<y<<"\t";
		}
		cout<<endl;
	}

	delete [] pivot;
	for(i=0;i<n;i++) delete [] a[i];
	delete [] a;

	return 0;
}
